import holoviews as hv
hv.extension('bokeh')
import numpy as np
import scipy.signal

2. Espectro y Transformada de Fourier

2.1. Espectro de frecuencia

2.1.1. Breve historia

  • Isaac Newton llamó espectros a los componentes que forman la luz blanca y que normalmente no se pueden ver

  • Newton mostró usando prismas que la luz blanca puede descomponerse en colores y viceverza

  • Hoy entendemos que la luz como onda tiene una frecuencia asociada y que cada color es una frecuencia particular

../../_images/fourier-newton.jpg ../../_images/fourier-prism.jpg

Paradójicamente, Newton nunca acepto que esto se debía a la frecuencia de la radiación ya que creía en la teoría corpuscular de la luz

2.1.2. Ondas y señales sinusoidales

Una onda es una perturbación que transporta energía a través del espacio

Tipicamente se describe por su frecuencia, amplitud y desfase. La frecuencia (Hertz) es el recíproco del período (segundos)

La siguiente es una función del tiempo completamente descrita por su amplitud \(A\), frecuencia \(f\) (período \(P=1/f\)) y fase \(\phi\)

\[ s(t) = A \cos (2 \pi f t + \phi) \]

Observe en el siguiente ejemplo como cambia la señal al modificar sus parámetros

def wave(time, frequency, phase, amplitude):
    return amplitude*np.cos(2.0*np.pi*time*frequency + phase)
signal_plot = hv.HoloMap(kdims=['Amplitud', 'Frecuencia', 'Desfase'])

t = np.linspace(0, 2, num=1000)
for A in [0.5, 1.]:
    for f in [0.5, 1, 2]:
        for phi in np.arange(0., np.pi, np.pi/10):
            signal_plot[(A, f, phi)] = hv.Curve((t, wave(t, f, phi, A)))
            
signal_plot.opts(hv.opts.Curve(width=500, height=300))

2.1.3. Componentes frecuenciales y descomposiciones armónicas

Una sinusoide que es períodica en una fracción

\[ \frac{P}{k} ~ \forall k \in \mathbb{N} \]

también es periódica en \(P\)

En general llamamos a

  • \(f_0 = 1/P\) la frecuencia fundamental

  • \(f_k = k/P = kf_0\) el k-esimo armónico de \(f_0\)

Si sumamos armónicos con distintas amplitudes el resultado es una nueva función periódica que tiene la misma frecuencia fundamental

En el siguiente ejemplo modifique las amplitudes de cada armónico y observe la señal resultante

\[ s(t) = A_1 \cos(2\pi f_0 t ) + A_2 \cos(2\pi 2 f_0 t) + A_3 \cos(2\pi 3 f_0 t) \]

Note que la suma de tres sinusoides no es necesariamente otra sinusoide

t = np.linspace(0, 2, num=1000)

def my_signal(time, frequency, Amplitudes):
    s = 0
    for k, A in enumerate(Amplitudes):
        s += A*np.cos(2.0*np.pi*time*(k+1)*frequency) 
    return s
signal_plot = hv.HoloMap(kdims=['Amplitud fundamental', 
                                'Amplitud 1er armónico', 
                                'Amplitud 2do armónico', 
                                'Frecuencia'])
for f in [1, 2]:
    for A1 in [0.5, 1.]:
        for A2 in [0, 0.5, 1.]:
            for A3 in [0, 0.5, 1.]:
                signal_plot[(A1, A2, A3, f)] = hv.Curve((t, my_signal(t, f, [A1, A2, A3])))
                
signal_plot.opts(hv.opts.Curve(width=500, height=300))

La noción de sintetizar señales en base a sinusoides se generaliza por medio de la serie trigonométrica

\[\begin{split} \begin{align} s(t) &= \sum_{k=0}^\infty A_k \cos(2\pi k f t + \phi_k) \nonumber \\ &= \sum_{k=0}^\infty a_k \cos(2\pi k f t) + b_k \sin(2\pi k f t), \nonumber \end{align} \end{split}\]

donde \(a_k = A_k \cos(\phi_k)\) y \(b_k = -A_k \sin(\phi_k)\) se obtienen de \(\cos(x+y) = \cos(x)\cos(y) - \sin(x)\sin(y)\)

Con esto podemos generar funciones periódicas arbitrarias definiendo \(\{a_k, b_k\}\) y \(f\)

Ejemplo

Sea por ejemplo \(a_k = 0\) y \(b_k = 1/k\) \(\forall k\)

¿Qué señal se obtiene al agregar cada vez más armónicos?

t = np.linspace(0, 5, num=1000)

def trigonometric_series(time, frequency, K):
    s = 0
    for k in range(1, K):
        s += np.sin(2.0*np.pi*time*k*frequency)/(k) 
    return s
signal_plot = hv.HoloMap(kdims=['Cantidad de armónicos', 'Frecuencia'])
for f in [1, 2]:
    for K in [1, 2, 3, 5, 10, 20, 50, 100]:     
        signal_plot[(K, f)] = hv.Curve((t, trigonometric_series(t, f, K)))
                
signal_plot.opts(hv.opts.Curve(width=500, height=300))

Nota

Hasta ahora hemos visto como sintetizar señales a partir de sus armónicos, pero también podemos hacer el proceso inverso, es decir encontrar los armónicos de una señal dada

Tal como Newton descompuso la luz blanca en colores, nosotros podemos descomponer una señal en sus armónicos usando la serie de Fourier

2.2. Serie de Fourier

Un poco de historia: En 1807 Jean Baptiste Joseph Fourier presenta un teorema indicando que una función periódica arbitraria con periódo \(P=1/f_0\) puede representarse como una suma ponderada de senos y cosenos

La serie de Fourier (FS) es una generalización de la serie trigonométrica a los números complejos. Puede revisar el apéndice de esta lección si necesita refrescar la memoria sobre los números complejos y sus representaciones

La FS se define como

\[ s(t) = \sum_{k=-\infty}^{\infty} c_k e^{j 2\pi k f_0 t}, ~~ c_k \in \mathbb{C} \]

A continuación veremos como obtener los coeficientes de Fourier \(c_k\) para una señal periódicas arbitraria

2.2.1. Base de Fourier

El conjunto de funciones

\[ v_k (t) = \frac{1}{\sqrt{P}} e^{j2\pi k t / P} ~~ \forall k \in \mathbb{Z} \]

se conoce como base de Fourier y cumple con

\[\begin{split} \langle v_n (t), v_m (t) \rangle = \int_0^P v_n (t) v_m^* (t) dt = \frac{1}{P} \int_0^P e^{j2\pi (n-m)t/P} dt =\begin{cases}1 & n=m \\ 0 & n \neq m\end{cases} \end{split}\]

La base de Fourier es un conjunto ortonormal en el espacio de funciones periódicas con periódo \(P\). Esto es facilmente comprobable si estudiamos

\[ \int_0^P e^{j2\pi k t / P} dt = \int_0^P \cos(2\pi k t/P) dt + j \int_0^P \sin(2\pi k t/P) dt \]

Podemos usar esta propiedad y escribir el producto punto entre una señal y el m-esimo elemento de la base de Fourier

\[\begin{split} \begin{align} \langle s(t), e^{j 2\pi m t/ P} \rangle &= \int_0^P s(t) e^{-j 2\pi m t/ P} dt \nonumber \\ &= \int_0^P \sum_{k=-\infty}^{\infty} c_k e^{j 2\pi k t/P} e^{-j 2\pi m t/P} dt \nonumber \\ &= \sum_{k=-\infty}^{\infty} c_k \int_0^P e^{j 2\pi (k-m) t/P} dt \nonumber \\ &= c_m P \nonumber \\ \end{align} \end{split}\]

Esto nos da una forma sencilla para obtener los \(c_m\) de una señal periódica arbitraria

\[ c_m = \frac{1}{P} \int_0^P s(t) e^{-j 2\pi m t/P} dt \]

2.2.2. Ejemplo formativo: FS de señal cuadrada

Sea

\[\begin{split} s(t) = \begin{cases} 1 & t \in[0, \frac{P}{2}] \\ 0 & t \in [\frac{P}{2}, P] \end{cases} \end{split}\]

Los coeficientes de su FS son

\[ c_0 = \frac{1}{P} \int_0^P s(t) dt = \frac{1}{P} \int_0^{P/2} dt = \frac{1}{2} \]

y

\[\begin{split} \begin{align} c_k &= \frac{1}{P} \int_0^\frac{P}{2} e^{-j2\pi kt/P} dt \nonumber \\ &= -\frac{j}{P} \int_0^\frac{P}{2} \sin(2\pi kt/P) dt \nonumber \\ &= 0 + j \frac{\cos(\pi k) - 1}{2\pi k} \nonumber \end{align} \end{split}\]

Notemos que a excepción de \(c_0\) los coeficientes sólo tienen parte imaginaria. Además sólo los armónicos impares (senos) son distintos de cero

Finalmente la FS está dada por

\[\begin{split} \begin{align} s(t) &= \sum_{k=-\infty}^{\infty} j \frac{\cos(\pi k) - 1}{2\pi k} e^{j 2\pi k t/P} \nonumber \\ &= \frac{1}{2} + \sum_{k=1}^{\infty} \frac{1 - \cos(\pi k)}{\pi k} \sin(2\pi k t/P) \nonumber \end{align} \end{split}\]

Programemos esta expresión y visualicemos el resultado al incrementar cada vez la cantidad de armónicos

t = np.linspace(0, 5, num=1000)

def FS_square(tiempo, frecuencia, K):
    s = 0.5 # C_0
    for k in range(1, K):
        c_k = (1-np.cos(np.pi*k))/(np.pi*k)
        s += c_k*np.sin(2.0*np.pi*tiempo*k*frecuencia)
    return s
signal_plot = hv.HoloMap(kdims=['Cantidad de armónicos', 'Frecuencia', ])

for f in [1, 2]:
    for K in [1, 2, 3, 5, 10, 20, 50, 100]:     
        signal_plot[(K, f)] = hv.Curve((t, FS_square(t, f, K)))
                
signal_plot.opts(hv.opts.Curve(width=500, height=300))

2.2.3. Propiedades de la FS

A continuación se enumeran algunas de las propiedades más importantes de la FS

  • Si \(s(t)\) es par entonces \(c_k\) es par

  • Si \(s(t)\) es impar entonces \(c_k\) es impar

  • Si \(s(t + P/2) = -s(t)\) (antiperiódica) entonces \(c_k=0\) para k par

  • Si \(s(t)\) es real y par entonces \(c_k\) es real y par

  • Si \(s(t)\) es real e impar entonces \(c_k\) es imaginario e impar

  • La FS es lineal

El Teorema de Parseval relaciona la potencia de una señal periódica arbitraria con sus coeficientes de Fourier

\[\begin{split} \begin{align} P_s &= \frac{1}{P} \int_0^P |s(t)|^2 dt \nonumber \\ &= \frac{1}{P} \int_0^P |\sum_{k=-\infty}^{\infty} c_k e^{j 2\pi k t/P }|^2 dt \nonumber \\ &= \frac{1}{P} \int_0^P \sum_{k=-\infty}^{\infty} |c_k |^2 dt \nonumber \\ &= \sum_{k=-\infty}^\infty |c_k|^2 \nonumber \end{align} \end{split}\]

2.3. Transformada de Fourier

Podemos extender el concepto de descomposición armónica a señales no periódicas usando la Transformada de Fourier. Esta herramienta matemática será fundamental en este curso

  • El concepto de frecuencia puede aplicarse también a señales no-periódicas

  • Según Joseph Fourier una señal no-periódica puede ser vista como una señal periódica con un período infinito

  • El único requisito es que ahora las frecuencias son un continuo, con un espaciado infinitesimal

A continuación veremos que una señal analógica puede ser vista como continua en el tiempo o continua en frecuencia. En la próxima lección veremos la Transformada de Fourier discreta (DFT) para señales digitales

2.3.1. Derivación de la transformada de Fourier a partir de la FS

Sea un tren de pulsos cuadrado con periódo P y ancho \(2T < P\) definido en un período como

\[\begin{split} s(t) = \begin{cases} 1, & |t| < T \\ 0, & T<|t| < P/2 \end{cases} \end{split}\]

Los coeficientes de su serie de Fourier son

\[ c_0 = \frac{1}{P} \int_{-P/2}^{P/2} s(t) dt = \frac{1}{P} \int_{-T}^{T} dt = \frac{2T}{P} \]

y

\[\begin{split} \begin{align} c_k &= \frac{1}{P} \int_{-T}^{T} e^{-j2\pi kt/P} dt \nonumber \\ &= \frac{1}{\pi k} \sin \left (2\pi k \frac{T}{P} \right) \nonumber \end{align} \end{split}\]

En el siguiente ejemplo la figura izquierda muestra la señal de pulso cuadrado mientras que la derecha muestra \(c_k\) en función de \(k\)

¿Qué ocurre con \(c_k\) a medida que \(P\) aumenta (\(f\) disminuye)?

def fs_coef(P, T, maxK=200):
    K = np.arange(-maxK, maxK)
    ck = np.zeros_like(K, dtype='float32')
    ck[len(K)//2] = 2*T/P
    for k in range(1, K[-1]):
        ck[len(K)//2+k] = np.sin(2*np.pi*k*T/P)/(np.pi*k)
        ck[len(K)//2-k] = ck[len(K)//2+k]
    return K, ck

def synthesis(time, ck, P):
    maxK = len(ck)//2
    s = ck[maxK]
    for k in range(1, maxK):    
        s += 2*np.cos(2*np.pi*time*k/P)*ck[maxK+k]
    return s   
t = np.linspace(-3, 3, num=2000)
T = 0.25

signal_plot = hv.HoloMap(kdims=['Periodo'])
fs_plot = hv.HoloMap(kdims=['Periodo'])
for P in [1, 2, 3, 5, 10]:
    K, ck = fs_coef(P, T)
    fs_plot[P] =  hv.Spikes((K, P*ck), 'Frecuencia (k)', 'P x ck').opts(xlim=(-50, 50))
    signal_plot[P] =  hv.Curve((t, synthesis(t, ck, P)), 'Tiempo (segundos)', 'Señal')
    

(signal_plot+fs_plot)

De la figura tenemos que cuando P es grande (o equivalentemente f es pequeño)

  • el tren de pulsos tiende a un único pulso

  • los coeficientes toman una “forma suave”

Llamaremos a esta forma suave de \(P \cdot c_k\) la “envolvente” de los coeficientes. Si tomamos el caso \(P \to \infty\) la envolvente resultante es

\[ \lim_{P\to \infty} P c_k = S(f) = \int_{-\infty}^{\infty} s(t) e^{-j 2\pi t f} dt, \]

o también

\[ S(\omega) = \int_{-\infty}^{\infty} s(t) e^{-j\omega t } dt, \]

donde \(\omega = 2\pi f\) se llama frecuencia angular.

Nota

Esto se conoce como transformada de Fourier directa o integral de Fourier

Reemplazando \(c_k = f_0 S(k f_0)\) y tomando el límite cuando \(f_0 = \frac{1}{P} \to 0\) tenemos

\[\begin{split} \begin{align} s(t) &= \lim_{f_0 \to 0} \sum_{k=-\infty}^{\infty} f_0 S(k f_0) e^{j 2\pi t k f_0} \nonumber \\ &= \int_{-\infty}^{\infty} S(f) e^{j 2\pi t f} df \nonumber \end{align} \end{split}\]

o también

\[ s(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} S(\omega) e^{j \omega t } d\omega, \]

que se conoce como la transformada de Fourier inversa

2.3.2. Par de Fourier

Los operadores (transformada de Fourier directa)

\[ S(\omega) = \mathbb{FT} [s(t)] = \int_{-\infty}^{\infty} s(t) e^{-j\omega t } dt \]

y (transformada de Fourier inversa)

\[ s(t) = \mathbb{FT}^{-1} [S(\omega)] = \frac{1}{2\pi} \int_{-\infty}^{\infty} S(\omega) e^{j \omega t } d\omega, \nonumber \]

se conocen como par de Fourier y nos permiten analizar una señal en el dominio del tiempo o en el dominio de la frecuencia sin pérdidas

2.3.3. Ejemplo: Transformada del pulso cuadrado

Consideremos nuevamente el pulso cuadrado

\[\begin{split} s(t) = \begin{cases} 1, & |t| < T \\ 0, & |t| > T\end{cases} \end{split}\]

su transformada de Fourier es

\[\begin{split} \begin{align} S(\omega) &= \int_{-\infty}^{\infty} s(t) e^{-j\omega t } dt \nonumber \\ &= \int_{-T}^{T} e^{-j\omega t } dt \nonumber \\ &= \frac{1}{-j\omega} \left(e^{-j\omega T } - e^{j\omega T } \right) \nonumber \\ &= \frac{2}{\omega} \sin(\omega T) = 2T \text{sinc}(\omega T) \nonumber \end{align} \end{split}\]

Esto calza precisamente con la “envolvente” que vimos anteriomente

T = 0.25

w = np.arange(-100, 100)
S = 2*np.sin(w*T)/w
S[w==0] = 0.5
<ipython-input-13-6be2a149569f>:4: RuntimeWarning: invalid value encountered in true_divide
  S = 2*np.sin(w*T)/w
hv.Curve((w, S), 'Frecuencia angular', 'Espectro').opts(width=500)

2.3.4. Propiedades de la transformada de Fourier

  1. La transformada de Fourier es un operador lineal, si tenemos dos señales y dos valores escalares entonces

\[ \mathbb{FT}[c_1 s_1(t) + c_2 s_2(t)] = c_1\mathbb{FT}[s_1(t)] + c_2\mathbb{FT}[s_2(t)] \]
  1. El operador de convolución en el tiempo se convierte en multiplicación en frecuencia

\[ \mathbb{FT}[(s_1 * s_2)(t)] = \mathbb{FT}[s_1(t)] \cdot \mathbb{FT}[s_2(t)], \]

donde

\[ (s_1 * s_2)(t) = \int s_1(\tau) s_2(t-\tau) d\tau \]

es la operación de convolución

  1. Así mismo, la multiplicación en frecuencia se convierte en multiplicación en el tiempo

\[ \mathbb{FT}[s_1(t)\cdot s_2(t)] = \frac{1}{2\pi}\mathbb{FT}[s_1(t)] * \mathbb{FT}[s_2(t)] \]

Finalmente el Teorema de Parseval nos dice que la energía se preserva entre frecuencia y tiempo

\[ \int | s(t) |^2 dt = \frac{1}{2\pi} \int | S(\omega) |^2 d\omega \]

2.3.5. Definición: Espectros de amplitud y fase

Llamamos a \(S(\omega)\) el espectro (transformada de Fourier) de una señal \(s(t)\)

El espectro es un número complejo que podemos escribir en notación polar como

\[ S(\omega) = |S(\omega)| e^{j\Phi(\omega)}, \]

donde \(|S(\omega)|\) se conoce como espectro de amplitud y \(\Phi(\omega)\) como espectro de fase

2.4. Resumen de la lección

En esta lección aprendimos

  • a realizar síntesis de señales usando sinusoides

  • a descomponer señales periódicas usando la serie de Fourier

  • a descomponer señales periódicas y no periódicas usando la transformada de Fourier

2.5. Apéndice: Números complejos

Sea z un número complejo, lo podemos escribir en forma cartesiana

\[ z = \Re[z] + j \Im[z] = a + j b \]

donde \(a \in \mathbb{R}\), \(b \in \mathbb{R}\) y \(j = \sqrt{-1}\) es el número imaginario.

También podemos escribirlo en forma polar

\[ z = c e^{j\phi} = c \cos(\phi) + j c \sin(\phi) \]

donde

  • \(c = |z| = \sqrt{a^2 + b^2} \in [0, \infty]\) es la magnitud

  • \(\phi = \angle z = \tan^{-1} \left (\frac{b}{a} \right) \in [-\frac{\pi}{2}, \frac{\pi}{2}]\) es el ángulo

  • \(a = c \cos(\phi)\)

  • \(b = c\sin(\phi)\)

se pueden escribir las siguientes relaciones

\[ \cos(\phi) = \frac{1}{2} (e^{j\phi} + e^{-j\phi}) ~\wedge~ \sin(\phi) = \frac{1}{2j} (e^{j\phi} - e^{-j\phi}) \]

el complejo conjugado de \(z = a + j b = c e^{j\phi}\) es

\[ z^* = a - jb = c e^{-j\phi} \]